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FERTUBBATION  SOLUIIONS  FOR  VARIABLE  ENERGY  BLAST  WAVES 


Edvard  T.  Pitkin 


ABSTRACT 

The  trajectory  of  and  the  flow  field  behind  blast  waves  with  tine  vaxying 
energy  input  is  determined.  Freeman's^  Lagrangean  coordinate  fomnilatlcn  is 
modified  to  Include  both  the  geometric  factor^  a,  for  plane,  cylindrical  and 
spherical  shocks  and  edso  non-integer  values  of  S,  the  energy  input  parameter, 
in  a single  conputationed  algorithm.  Numerical  problems  associated  with  vanish- 
ing density  at  the  fictitious  piston  face  are  then  examined  and  solved.  Second 
order  perturbation  solutions  about  the  infinite  strength  shock  are  then  ob- 
tained in  Sakural's^  inverse  shock  Mach  number  expansion  parameter  for  0 S ^ < 
a 4-  1.  Tables  and  graphs  of  sl0:ilfic8nt  numerlce^.  coefficients  are  presented 
for  cawpariBaa  to  and  extension  of  results  of  other  authors.  Graphs  of  Igrplceil 
shock  trajectories  and  flow  field  density,  pressure  and  velocity  variations 
are  also  presented  and  discussed. 
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1. 

PERTURBATION  SOLUTIONS  FOR  VARIABLE  ENERGY  BLAST  WAVES 
...  t 

Edward  T.  Pitkint 

• • / 

1.  INTRODUCTION 

The  analytical  theory  of  blast  waves  hats  been  the  subject  of  intensive 

a 

study  by  many  investigators  since  the  pioneering  work  of  G.  I.  Taylor  about 
thirty  years  ago.  A large  fraction  of  this  work  concerns  waves  of  constant, 
instantaneously  supplied  energy  which  propagate  and  decay  according  to  the 
dictates  of  inviscid  isentropic  flow  behind  the  wave  and  Rankine-Hugoniot 
conditions  at  the  wave  front.  This  work  has  been  well  summarized  by  Sakurai^ 
who  has  also  developed  an  efficient  method  of  expanding  the  general  solution 
filjout  ttio  olmpLer  limiting  solution  for  the  flow  behind  a very  strong  shock. 
Furl.hcrmore , he  has  shown  that  the  square  of  the  inverse  shock  Mach  nvimber  is 
a most  convenient  and  natural  parameter  to  use  for  this  expansion. 

A lesser  amoimt  of  work  has  been  devoted  to  analysis  of  blast  waves  of 
variable  energy  input.  In  1957  Lees  and  Kubota®  examined  this  case  briefly 
in  connection  with  the  hypersonic  blast  wavs  analogy  and  Rogers^  presented 
the  limiting  strong  shock  solution  in  1958.  Mcy-;  .'r  '-iutly  in  1968  Freeman^ 
examined  this  case  in  connection  with  cylindrical  spa.^-k  channel  formation 
from  exploding  wires  and  Dabora^  has  applied  variable  blast  wave  techniques 
to  droplet  explosions  in  spray  detonations. 

'Phi!  work  reported  hore  is  bauicully  un  extension  of  Freeman's  approach 
to  cover  the  full  spectrum  of  plane,  iTylincirical  and  spherical  variable 

t Professor  of  Mechanical  and  Aerospace  Engineering,  The  University  of 
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energy  blast  waves  within  a single  analytic  formulation  and  computation  algo- 
rithm. It  also  contains  corrections  to  some  numerical  errors  discovered  in 
Freeman's  paper.  As  a consequence  Freeman's  notation  has  been  adopted  com- 
pletely and  derivations  given  in  his  paper  have  only  been  summarized  here. 

His  major  contribution  was  to  recast  Sakurai's  expansion  procedtire  in  La- 
grangean  coordinates,  replacing  distance  from  the  origin,  r,  with  mass,  m,  as 
an  independent  vsiriable,  but  retaining  the  inverse  shock  Mach  number  as  the 
expansion  psirameter.  The  advantage  is  two  fold;  (l)  the  number  of  differ- 
ential eq\iations  that  must  ultimately  be  integrated  numerically  is  reduced  by 
nearly  a factor  of  two,  and  (2)  the  position  of  the  fictitious  m'^ving  piston  face 
which  supplies  the  variable  energy  is  uniquely  defined  to  be  at  mass  = 0.  In 
Kulerian  coordinates  precise  determination  of  the  piston  position  can  be  a 
very  difficult  numerical  task. 

In  sections  2 and  3,  the  equations  of  motion  emd  the  strong  shock  solution 
are  reviewed.  The  expansion  procedure  is  outlined  in  section  i*  and  a method 
of  extrapolating  the  numerical  solution  to  the  mass  origin  is  discussed  in 
section  5.  Numericai  results  are  then  summarized  and  discussed  in  section  6. 


2.  EQUATIONS  OF  MOTION 

I 

Let  Pq  and  Pq  be  the  pressvire  and  density  in  the  vindisturbed  medium 
while  r is  the  radius  from  the  origin  to  an  arbitrary  point  in  space  and  M is 

the  total  mass  contained  between  the  origin  and  r.  Let  r = R at  the  shock 

front  so  that  the  volume  behind,  the  shock  front,  when  filled  with  gas  of  orig- 
inal density,  will  contain  the  mass 

M„  = ' (2.1) 

where  a is  a geometric  parameter  with  values  0,  1,  and  2 for  plane,  cylindri- 
cal and  spherical  waves  respectively.  Finally  let  u be  the  speed  at  any 

point  and  U be  the  speed  of  the  shock  front,  C the  speed  of  sound  and  Y the 

specific  heat  ratio  in  the  undistiirbed  gas. 

It  is  convenient  to  introduce  non-dimensional  variables.  For  mass,  rad- 


ius, velocity  and  density  teike 

ra  = M/Mo  (2.2) 

a ^ r/R  (2.3) 

f = u/U  (2.1+) 

h = p/Po  (2.5) 

while  the  sqiiare  of  the  inverse  shock  Mach  number 

y = C2/U2  = c2/(dR/dt)2  (2.6) 


may  be  combined  with  to  define  a variable  related  to  non-dimensional  pres- 
:**urc 


g = py/Po 


(2.7) 


An  arbitrary  non-zero  initial  radius  of  the  shock,  R^,  is  then  used  to  define 


a characteristic  time 


^ (2.8) 

which. may  then  be  used  to  non-dimensionalize  time  through 

0)  = t/to  (2.9) 

while  the  shock  front  coordinate,  R,  is  non-dimensionalized  by 

z = R/Ro  (2.10) 

Finally  the  decay  index,  a ratio  of  the  percentage  changes  of  shock  Mach  num- 
ber and  radius,  as  introduced  by  Sakurai^,  is 


X = = 4R£SZd^ 

Tdz/z)  (dR/dt)2 


(2.11) 


the  latter  obtained  by  use  of  (2.6)  and  (2.10). 

The  equations  of  fluid  flow  behind  the  shock  wave  have  been  derived  in 
Lagrangean  coordinates  by  Freeman^.  Only  the  results  will  therefore  be  sum- 
marized here.  First  the  density  at  any  point  is  given  by 

h = p/Po  = l/[(a+l)a'^3xi/3m]  (2.12) 

and  the  non-dimensional  velocity  is 

f = a-(a+l)m|~  + ^y^'*  (2.13) 

Beginning  with  the  assumption  of  a pxurely  isentropic  process  behind  the  shock 
for  each  mass  element  it  follows  that 

8(p/p^)/8y  = 0 (2.1U) 


- U - 


' r 


which  upon  iiaing  the  non-dimenalonal  variables ^ above  can  be  shown  to  lead  to 
the  following  in  Lagrangean  coordinates , 

(«.♦!)«  ria  |4  t i Ja)  M t » i M 

L\d  3m  Yg  • 9®  3m2J  y 3m 


, ria  3d  ^ 1 ^ 3^d  1 ^ - 

^ lid  3y  Yg  3y ' 3m  3m3yJ  • 


(2.15) 


Similarly  Newton's  second  law  for  the  mass  element 


3u  _ 3u  dv 
3t  3y  dt 


-(2r)V(3-“^/23| 


loud::  to 


/l“(a+l)  3g  Xd  . / . , s2o+X  id  . / ..\2  oi^d 

is  - T * * <“*!> 

, r O/  \ -.3^d  . 2+X  3d  . 3a  3X  . . 3^d1  . 


(2.16) 


(2. IT) 


The  energy  contained  between  the  origin  and  the  shock  is  the  sum  of  the  kine- 
tic energy  and  the  internal  energy.  The  increment  of  energy  over  that  of  the 
same  volume  of  undisturbed  gas  is 


^Mo 

■n 


Po 


^ ^ 2*^  (y-i)p  ~ (y-1)p, 


•|  dM  , 


(2.18) 


Converting  to  non-dimensional  form,  using  Equation  (2.12)  and  letting  a prime 
designate  D/3m  gives 
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(a+l)p 


1 r / Yf^  . \ 


(2.19) 


Letting  b = 1/[(y-1)(oi+1)  ]*  j“  and  designating  the  integral 

J after  Sakurai^  this  becomes 


Eg  _ J , 


(2.20) 


If  the  total  energy  behind  the  shock  varies  as  a pwer  of  time  according  to 


the  following  rule 


^g  = Ro“^^Po(t/to)^ 


the  energy  equation  reduces  to  the  simple  form 


f?  / (1+1  T / -u 

w /z  * J/y  - b. 


(2.21) 


(2.22) 


When  6=0  the  energy  is  constant,  a case  examined  by  many  investigators  and 
we’’!  documented  by  Sakurai.  Letting  6 take  non-zero  values  allows  one  to  model 
a fairly  wide  reinge  of  energy  inputs  which  can  be  further  extended  with  some 
modification  of  (2.2l)  as  has  been  shown  by  Freeman^. 
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3.  STRONG  SHOCK  LIMIT 

At  very  large  shock.  Mach  numbers  y 0 and  J/y  >>  b so  the  energy  equation 
reduces  to  a simple  form  consistent  with  Rogers'^  similarity  solution  for  the 


equation  of  motion  for  which  the  shock  front  advances  according  to 


R = = Kt" 


(3.1) 


where  n is  related  to  R,  a and  A through 


n = (2+6)/(a+3)  = 2/(Ao+2) 


(3.2) 


from  which  the  decay  index  for  strong  shocks  is 


A^  = 2(a+l-6)/(2+6) 


(3.3) 


The  zero  subscripts  on  A and  J are  used  to  designate  the  strong  shock  case 


here. 


Note  also  that 


^ - ill.  ^ U _ 1 
du  ■ Rg  dt  C " ^ 


(3.14) 


while  from  (3.l) 


dz/d(o  = nz/(o 


(3.5) 


so  that 


z = lo/ny* 


(3.6) 


Combining  this  with  the  strong  shock  energy  equation  (b  = 0)  yields 


(3.7) 
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4.  OUTLINE  OF  SERIES  SOLUTION 


Eqmtions  (2.15)  (2.17)  and  (2.22)  comprise  a system  of  three  dependent 
variables;  A,  a,  and  g and  two  independent  variables;  m and  y.  Sakurai  has 
shown  that  A is  a function  of  y only  and  that  the  pressure,  density  etc.  are 
slowly  varying  functions  of  y.  He  therefore  suggests  expsuasion  in  powers  of 
y to  reduce  the  number  of  independent  variables.  In  the  Lagrangean  coordinates 
used  here  this  has  the  very  desirable  effect  of  reducing  the  partial  differ- 
ential equations  to  ordinary  differential  equations. 

/ 

The  expansions,  taken  about  the  strong  shock  solution  for  which  y = 0, 

are 


A(y)  = 

A + A y + A y^  + . . . . 

(4.1) 

o 1 2 

d(m,y) 

= aQ(m)  + aj(m)y  + a2(“)y‘^  + . . . . 

(4.2) 

g(ni,y) 

= RqM  + gj(m)y  + g2(ni)y^  + . . . . 

(4.3) 

By  similar  reasoning  the  equation  for  the  shock  front  is  expanded  about  the 
strong  shock  "zero  order"  solution  (3.7)  as 

z = Zq  (l  + Zjy+Z2y^  + ***)  • (4.4) 

Equation  (3.4)  which  applies  in  general  can  be  integrated  to  give: 


Substituting  for  z and  taking  z = 0 at  u = 0 then  gives  an  expansion  for  the 
time  variable  as  a function  of  y and  z 


u = 


+ “jZjY  + •••  + UjZjyJ  + •••) 


(U.6) 


where 


Uj  :=  (2jAQ+2)/((2J+l)XQ+2)  / 


(U.7) 


Noting  from  (2.13)  that  f = f(A,  a,  m)  it  is  clear  that  substitution  of 
(U.l  - U.3)  into  (2.19)  will  yield  integrals  in  ascending  powers  of  y so  that 
J can  be  expressed  as 


J = Jo 


(1+.8) 


These  expansions  are  then  substituted  into  the  isentropic,  momentum  and 
energy  equations  (2.15),  (2.1?)  and  (2.22)  and  the  coefficients  of  like  powers 
of  y equated  to  give  perturbation  equations  of  ascending  order.  In  the  isen- 
tropic and  momentum  equations  this  procedure  yields  ordinary  differential 
equations  for  aj  and  gj  with  Aj  entering  as  a parameter.  These  must  be  solved 
subject  to  the  shock  Jump  boundary  conditions. 


dod)  = 1 


g^d)  = 2^/(y+1) 


djd)  . 0 


<(i) 


Y - 1 

(a+l)(Y+l) 


gj(l)  = (i-y)/(y+i) 


■ (a+l)(Y+l) 


(d9) 


The  remaining  coefficients  are  all  zero  at  the  shock.  Fortimately  the 
successive  orders  of  approximation  to  the  isentropic  equation  yield  intermedi- 
ate integrals  for  gj.  These  may  then  be  used  in  the  momentum  eqiiation  which 
must  be  solved  by  numerical  integration  from  the  shock  front  to  the  mass 
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origin.  The  equations  and  solutions  for  each  order  will  be  discussed  in  the 
following  sections. 

Now  consider  the  energy  equation.  Putting  the  expansion  (U.U),  {U.6) 
and  {U.8)  into  (2.22)  gives 


= Zo^^(l  + ZjY  + '''  ***^' 


The  two  series  with  exponents  may  be  converted  to  simple  power  series  by  means 
of  the  binomled  formula  if  the  first  terra  of  each  is  larger  than  the  sum  of  the 
remaining  terms.  This  is  certainly  true  for  small  y.  Multiplying  out  on  the 
right  and  equating  coefficients  of  like  powers  of  y then  yields  relations  be- 
tween Jj  and  Zj  ; 


+ rii£2.  . (a+i)]z,  (U.ll) 

Luq  J 2 • 


These  relations  are  eqxially  valid  for  integer  and  non  integer  values  of  6. 
Freeman^  only  considered  integer  values.  Substituting  (U.U)  into  (2.11)  and 
solving  for  the  Zj  gives 


r- 


^1  ~ 

Z,  = -[X,  + (l+X^)z,X  l/2xj 
^ ^ olio 

• i 

Z3  = -[X3  + il+X^)z^X2  + (1+2Xq)z2X^]/3a2  (U,12) 

» • 

# • 

* ♦ 

which  may  then  he  used  in  (U.ll)  to  give  relations  between  the  Xj  and  Jj. 

These  will  be  used  later  to  evaluate  the  Xj  for  each  order  of  approximation. 

The  solution  proceeds  in  ascending  orders  of  approximation.  The  zeroth 
order,  corresponding  to  the  strong  shock  solution  with  y 0,  is  used  in  the 
computation  of  the  first  order  solution  and  the  first  order  solution  is  then 
used  1.0  compute  the  second  and  so  forth.  In  this  paper  the  solution  will  be 
curried  only  to  the  2nd  order. 

The  zeroth  order  approximation  to  the  isentropic  equation  is 


(it. 13) 


which  Freeman  has  shown  to  have  an  integral 


r,  - V /in  “/TMYn,^o/(“+l) 

60  ^O'  ^ 


(U.lU) 


where 


" Y+l  [(a+l)(Y+l)] 


(i^.l5) 


The  momentum  equation  in  the  zeroth  approximation  becomes 


I 


J 
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whf.'re 


Aq  = (a+l)Y  l%&o/‘^o  ~ t 


= (ci+l)mY(a+X^/2)  , 


Co  = Y[Xq/2  + a(a+l)g^a“a^/a2]  . 


(1*.17) 


The  first  approximation  to  the  energy  integral,  J,  of  (U.8)  will  he  obtained  hy 
evaluating 


= 


X (ao-(tt-*-l)man)^  ^ So^nflo 


a+1 


¥] 


dm 


(i4.l8) 


along  with  Integration  of  (I+.I6)  from  the  shock  front  to  the  mass  origin. 

(Note  that  Freeman's  &rras  of  this  and  higher  order  integrals  are  only  correct 
for  0=1.)  The  first  term  of  this  expression  is  a measure  of  the  kinetic  en- 
ergy behind  the  shock  and  the  second  term  is  proportional  to  the  internal  ener- 
gy so  these  two  factors  can  easily  be  accmulated  separately.  The  integration 
of  aQ  is  initiated  with  the  boundary  conditions  (4.9).  The  values  of  aQ(m), 
^(m),  go("')  s^nd  go(m)  may  be  stored  for  use  in  the  first  order  solution  or, 

nore  conveniently^ these  integrations  can  be  performed  simultaneously  with  those 
for  the  higher  order  solutions. 

In  the  zeroth  order  equation  \q  appears  as  a known  parameter  obtained  from 
the  similarity  solution.  In  the  higher  order  equations  the  ^j's  appear  as  \in- 
known  parameters.  Fortunately,  the  equations  are  linear  in  Xj  and  of  course 
the  perturbation  equations  are  linear  differential  equations  so  it  follows  that 
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superposition  of  solutions  is  valid.  This  property  vrill  allow  evaluation  of 
the  Ij's  after  the  equations  have  been  integrated. 

In  each  higher  order  approximation  the  gj  term  in  the  momentum  equation 
may  be  eliminated  with  the  isentropic  equation  and  its  integral,  gj.  The  inte- 
grals havd  been  obtained  through  third  order  by  Preemsm.  The  momentum  equation 
then  gives  higher  order  perturbation  equations  of  the  form 


(4.19) 


in  which  the  coefficients  Bj  through  Ej  are  functions  only  of  known  quantities 
from  lower  order  solutions.  It  follows  that  can  be  expressed  as  the  linear 


combination 


- ^.11  * ^J^j2 


(4.20) 


so  that  (4.19)  gives  two  equations 


“Jl  = "jl  * “Jl  ♦ , 


“J2  = <“j  “52  * “J2  * • 


(k.si) 


It  is  convenient  to  let  the  non-zero  botindary  condition  of  (4.2)  be  satisfied 
by  so  the  remaining  perturbations  are  zero  at  the  shock.  Note  that  the 
denominator,  in  the  same  for  all  orders  of  approximation.  This  term  goes 
to  zero  at  m = 0 when  6=0  causing  computational  difficulties  which  will  be 
disBcussed  later. 

It  will  be  evident  upon  substitution  of  Ctj  into  the  isentropic  equation 
that  gj  must  be  of  the  form 


^ «J2  ' 


(4.22) 
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Finally  substituting  and  gj  into  the  energy  integral  will  show  Jj  to  be  of 
the  form 


•'j  ■ •'ji  * '/ja 


)c 


{‘♦.23) 


where  the  terms  Jj,  and  Jj2  are  to  be  accumulated  simultaneously  with  integra- 
tion of  ( ‘♦.21).  When  Jj^  and  Jj2  have  thus  been  evaluated  Xj  can  be  deter- 
mined. After  using  (‘♦.12)  in  (‘♦.ll)  to  eliminate  Zj  in  favor  of  Xj,  (‘♦.23)  is 

/ 

substituted  in  to  the  left  hand  side  of  (‘♦.ll)  and  the  result  is  then  solved 
for  X j . The  first  two  of  these  coefficients  are 


b - J 


X,  = 


11 


‘ J12  + Jo[a)j6/a)Q-(a+l)  1/X2 


(‘♦.21*) 


and 


X = Jq  M G - 2^0  J21 
Jq  II  - :>x2  J22 


(‘♦.25) 


where 


G = 6(6-1)  + (a+l)(a+l-Xo-26a)i/a)jj) 


+ (I+Xq)  6 0)2/01)0 


H = 6 (0)2/000) 


(‘♦.26) 


The  (i)j  are  obtained  from  (‘♦.7).  Expressions  for  the  coefficients  Aj  through 
Ej,  Jj^  and  Jj^  are  listed  in  the  Appendix. 

Freeman  has  shown  that  given  approximations  through  the  k'th  order  with 
k > 2,  it  is  possible  to  obtain  an  excellent  approximation  to  the  shock  posi- 
tion R(t)  as  follows.  Truncate  the  X series  at  the  k+2  term  and  evaluate  the 
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1 1.  1,1.1  ivii  11 


1^ 


laat  two  Aj  in  terms  of  those  preceeding  and  known  boundaries  of  the  X vs  y 
curve  at  y = 1.  Thus  ' 


X(l)  = Xq  + Xj  + •••  + Xk  + \+i  + ^k+2  ' 0 

I 

} 


(U.27) 


dx(l)/dy  = Xj  + 2X2  + •••  + kXj^  + (k+l)Xj^^^  + (k+2)Xj^^^ 

= (2+a)A  ('^•28) 

The  higher  order  X terms  in  (It.  12)  are  then  zero  and  all  the  Zj  cSp  be  expres- 
sed in  terms  of  X^  through  ^^+2*  If  k = 2 the  general  expression  for  Zj  when 
J 1 5 is 


2.  = -^[l+(j-lt)Xo]X^Zj_^  + [l+(j-3)X^]X,z 


“J 


+ [1+(J-2)X^]X2Z_j_2  + tl+(.)-l)X^]XjZj_J/jXo  (It. 29) 

Using  this  and  (H.  l?)  In  (li.lt)  and  (It.f))  then  gives  z(u))  or  R(t)  for  all  vaJuer. 


-if  y. 

It  was  found  that  this  procedure  works  very  well  for  3=0,  but  when  B 
approaches  a + 1 where  X = 0,  satisfaction  of  the  final  slope  condition  (It. 28) 
with  a limited  nximber  of  terms  causes  a false  hump  in  the  X-y  curve  near  y = 1. 
The  procedure  was  therefore  modified  through  multiplying  ^^+2  obtained  by 
simultaneous  solution  of  (It. 27)  and  (lt.28)  by  the  factor  [l-3/(a+l)]  and  solv- 
ing (It. 27)  for  X once  again.  Thus  the  need  to  meet  the  final  slope  condi- 
K. ' * 

tioris  Is  progressively  eliminated  as  3 approaches  a + 1 to  give  more  realis- 
I. in  A-y  curves. 


5.  EXTRAPOLATION  TO  ra  = 0. 


III  II, 


When  li  = 0 as  it  does  in  the  constant  energy  blast  vuv«,  integration  of 
(l*.l6)  will  show  that  ■*  os  and  <1^  ->•  0 as  m •>  0 so  it  follows  from  ('..17) 
that  -►  0.  This  latter  term  appears  as  the  denominator  of  {h.2l)  so  it  may 
be  anticipated  that  computational  difficulties  will  be  encountered  near  the 
mass  origin  for  all  orders  of  approximation.  The  problem  is  fxirther  compounded 

i 

in  evaluation  of  the  energy  integrals  which  contain  terms  proportional  to  a.^ 
and  accimulate  rapidly  near  m = 0.  For  example  if  Y - 1.1  and  a = 1,  sixty 
percent  of  the  true  value  of  Jq  is  accximulated  in  the  last  one  percent  of  mass 
variation.  NumerienL  integration  all  the  way  to  the  origin  is  clearly  impos- 
r.ihle  yet  it,  i r.  mandatory  that  accurate  values  of  Jj  be  obtained  else  the  next 
higiier  order  apiiroxiraation  be  seriously  in  error. 

Freeman^  mentioned  this  problem  and  suggested  an  "appropriate  approxima- 
tion formula"  but  gave  no  details.  Similar  numerical  problems  arising  in  the 
Eulerian  formation  of  the  problem  have  been  noted  by  Sakurai^  and  Bach  and 
Lee**.  The  procedure  outlined  here  was  developed  by  the  author  and  is  based 
upon  the  fact  that  in  a region  close  to  the  mass  origin  the  pressure  term,  g^, 
c;in  be  taken  to  be  constant  to  a very  high  degree  of  approximation.  Although 
tfie  procedure  if.  only  ne;cessary  when  3 is  close  to  zero,  it  is  still  valid  for 
larger  values  of  3 scj  it  con  be  incorporated  into  a computer  program  valid  for 


a I I values  of  3. 

Consider  the  integrand  of  (H.l8)  in  which  the  first  term  is  related  to 
the  kinetic  energy  and  the  second  to  the  internal  energy. 


2 ggaggp 


(5.1) 


,-3 


In  the  region  near  the  origin,  say  ra  < 10  , the  density  is  very  low  so 

it  may  be  anticipated  - and  verified  by  numerical  experiment  - that  the  first 
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term  is  much,  smaller  than  the  second.  Furthermore  it  has  been  observed  from 
numerical  integrations  that  the  product  ioAq  approaches  zero  at  nearly  the  same 
rate  as  so  that  this  first  term  is  nearly  constant.  These  numerical  inte- 
grations also  showed  g^  to  be  constant  to  about  four  significant  figure;;  in 
this  region.  It  follows  that  taking  these  factors  constant  over  an  interval  of 
dm  < 10  will  introduce  no  significant  error  in  the  conqputation  of  J^. 

The  major  contributing  factor  to  the  integral  is  then  which  by  use 

of  (lt.ll+)  can  be  expressed  in  terms  of  m as 


-(V  /rr  -XoA(«+l) 

^0^0  '^0'®0  ® 


(5.2) 


which  can  bo  inbegrated  analytically  if  g^  is  taken  to  be  constant. 

The'  procedure  than  1:;  to  lntf?grul.e  numericaJ  ly  from  rn  - 1 to  m = <S  where 
0 - iS  10  * and  add  a correction  term  obtained  by  analytical  Integration  J’rom 
0 to  6^  i .e. , 

1^^^  r;cl-^nA(a+l)' 


^ 2(a+l)^^o-(“‘*’l)^«o^^  l^j  1_Xo/y(o+1) 


(5.3) 


where  and  are  evaluated  at  m = 6. 


Similar  considerations  apply  to  the  higher  order  approximations  to  J. 

In  each  case  the  kinetic  energy  term  in  the  integrand  can  be  taken  as  constant 
over  the  last  integration  interval  while  the  internal  energy  term  can  be  ex- 
pressed in  terms  of  Bq  an*!  essentially  constant  ratios  such  us  a^/a^, 

a’]/aQ,  etc.  This  method  has  been  checked  by  letting  6 take  on  values  from 
lO”^  to  lO”^  and  noting  that  all  give  the  same  final  values  for  each  Jj . A 
value  of  6 = lO”^  was  used  in  the  computations  leading  to  the  results  reported 
in  the  next  section.  Formulas  for  the  correction  term  AJ^^,  etc.  are 

given  in  the  appendix. 
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6.  RESULTS 


The  eqmtions  outlined  in  Sections  3,  ^ and  the  Appendix  have  been  numeri- 
cally integrated  with  a fourth  order  Runge-Kutta  algorithm  in  double  precision 
and  extrax^olated  to  m = 0 as  discussed  in  Section  3.  The  procedure  has  been 
carried  out  to  the  second  order  in  the  expansion  variable  y (some  authors^ 
prefer  "third  approximation")  for  plane,  cylindrical  and  spherical  blast  waves, 
i.e.,  for  a = 0,  1,  2.  The  range  0 < 3 < “■*■!  (above  which  the  zeroth  order 

similarity  solution  ceases  to  exist^)  was  covered  for  Y - 1.1  and  I.I4.  The 
solution  for  the  case,  3 = a + 1,  which  corresponds  to  the  well  known  self  sim- 
ilar solution  for  which  A = 0 £ind  R = Kt,  is  computationally  unobtainable  be- 
cause of  a zero  divisor  problems  evident  in  equation  3.7.  However,  solutions 
for  3 very  close  to  this  value  (within  .01)  were  obtained  easily.  Solutions 
have  also  been  obtained  for  other  typical  values  of  Y for  integer  values  of  3. 

Tables  1 to  3 are  summaries  of  the  expansion  coefficients  Aj,  Xy  and 
Jo>  *^2  obtained  in  the  computation.  Korobienikov  and  Chushkin®,  Sakurai  * , 
Bach  and  Lee**  have  all  published  values  of  Aj  computed  in  Eulerian  coordinates 
for  3 = 0,  Y = l.U  which  agree  with  those  given  here.  (Note  that  Bach  and  Lee 
report  values  of  A^/S  while  the  others  give  values  of  Ax/(a+l)).  The  values  of 
A2  smd  A^  for  3=0  and  Y = l.U  are  also  in  agreement  to  the  full  six  signifi- 
cant figures  quoted  by  Bach  and  Lee.  The  values  of  A 1 for  3 = 0 and  Y = 1.1  to 
3.0  are  in  agreement  with  those  reported  by  Freeman,  however  his  values  of  Aj 
for  3=1  are  in  error  because  he  apparently  used  an  incorrect  value  of  Aq  to 
obtain  them.  His  values  of  Jjj  and  J12  are,  however,  correct.  As  s final  check 
on  numerical  accxiracy  it  is  noted  that  the  values  of  Jq  reported  to  1+  signifi- 
cant figxu-es  by  Rogers  for  Y = 1.2  and  1.4  and  various  3 agree  with  values  in 
these  tables. 

It  can  be  seen  from  the  tables  that  the  A coefficients  approach  zero 


rapidly  as  3 is  increased  from  zero.  This  is  also  shown  in  Figure  1.  The  max- 
imum value  of  3 is  a+1  where  X = 0 and  the  shock  propagates  with  constant  ve- 
locity. This  drop  off  of  the  decay  index  implies  much  stronger  convergence  of 
the  perturbation  series  as  3 increases. 

Figure  2 shows  typical  variations  of  the  J coefficients  with  3 for  « = 2, 
and  Y = 1.^4.  It  is  seen  that  J2  drops  off  rapidly  and  becones  fairly  con- 
stant while  Ji  takes  on  an  increasing  portion  of  the  energy  as  3 increases.  j 

Similar  results  obtain  for  other  Y and  a.  ' 

Radial  pressure,  density  and  velocity  destributions  are  presented  in 
Figures  3 to  5 for  u = 2,  Y = 1.1(  and  3 = 0,  1,  2.  The  hump  in  the  density 
distubution  for  y = .5  and  3 = 0 is  most  likely  due  to  truncation  of  the  expan- 
sion at  the  third  term.  This  suggests  an  upper  limit  of  validity  for  y = .5,  ' 


or  a shock  Mach  number  of  1.^,  for  the  second  order  expression  when  3=0. 

ITiis  limit  can  be  raised  considerably  as  3 increases  as  in  evidenced  by  the 
curves  for  3 = 1,  2. 

As  3 is  increased  from  zero,  the  constant  energy  case,  the  flow  region 
rapidly  contracts  to  a thin  shell  between  the  shock  and  the  fictitious  piston 
located  at  the  zero  mass  position  as  is  evident  in  Figures  I4  and  5.  Even  when 
3=0,  over  99%  of  the  mass  is  located  in  the  outer  half  radius  of  the  sphere 
but  as  noted  before,  over  30%  of  the  total  energy  can  be  confined  to  the  inner 
ifc  of  the  mass  because  the  temperature  becomes  extremely  large  in  this  region. 
As  3 approaches  a+1  the  solution  approaches  the  self  similar  result  for  X = 0 
and  y = constant.  At  this  limit  the  spherical  blast  wave  flow  field  occupies 
a thin  shell  within  which  the  velocity,  pressure  and  density  are  constant. 
Similar  results  obtain  for  other  Y and  a.  In  particular  the  profiles  for 
a = 1,  3 = 1 and  Y = 1.^  are  given  in  Figure  6 as  replacements  for  similar 


plots  given  by  Freeman  which,  as  noted  before,  are  in  error  due  to  his  use  of 
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the  incorrect  value  of  A. 

The  VEiriation  of  the  decay  index.  A,  and  the  shock  position,  z,  with  shock 
strength,  y,  is  given  in  Figure  7 for  a spherical  shock  with  Y = l.U  and  6=0. 

The  decay  index  was  obtained  by  use  of  accurate  values  of  A^,  Aj  and  A2  given 
in  Table  2 and  approximate  values  of  A3  and  A4  obtained  by  matching  boiuidary 
conditions  at  y = 1 as  outlined  in  section  U.  The  shock  position  was  obtained 
from  (U.6)  with  zj  from  (lt.l2)  and  enough  terms  of  (U.28)  to  insure  that  no 
truncation  error  would  occur  in  the  sixth  significant  figure.  These  results 
are  compared  with  a decay  index  and  shock  position  obtained  from  tables  of 
shock  position  and  over-pressure  reported  by  Goldstine  and  Von  Neumann®  who 
solved  the  exact  partial  differential  equations  by  a numerical  technique.  The 
comparison  is  most  favorable  over  the  whole  range  of  y even  though  the  internal 
structure  of  the  flow  field  can  only  bt;  predicted  to  y = .5  for  6=0.  Of 
(.•(jurse  when  6 becomes  larger  the  Internal  approximations  improve  for  larger  y. 

The  effect  of  6 on  decay  index  and  shock  position  is  given  in  Figures  8 
and  9 where  A and  z are  plotted  against  shock  strength  with  6 as  a parameter. 

Note  the  smooth  approach  to  A = 0 as  6 approaches  a+1  in  Figure  8. 

It  is  evident  in  Figure  9 that  a constant  value  of  shock  strength  is  ap- 
proached in  the  limit  6 = a+1  or  6 = 3 in  this  case.  The  value  can  easily  be 
obtained  by  noting  that  a solution  of  the  form  z is  required  in  which 

n = (2+6)/(a+3)  = 1 to  be  consistent  with  the  rest  of  the  family  of  ciurves. 

Differentiating  this  relation  and  using  (3.^),  the  definition  of  y,  gives 

d’/./du)  = z/u  = (6.i) 

ij 

while  the  energy  equation  (2.22)  gives 

(w/z)^^^  = Jy  - ]^[ (a+l) (y-I) ] (6.2) 


_ oi  _ 


liiiminating  u/z  and  using  the  second  order  approximation  for  J consistent  with 
the  other  curves  then  yields  the  following  equation  which  can  be  solved  for 
the  limiting  value  of  y. 


y(a+3)/2  = _ i/I(a+i)(Y-l)jj  y + J2  y^  (6.3) 

Taking  V = 1.*+,  the  limits  are  y = 0.531,  0.528  and  0.5^2  for  a = 2,  1 and  0, 

I 

respectively. 

The  shock  trajectories  corresponding  to  Figtires  8 and  9 are  given  in  Fig- 
ure 10.  It  is  clear  that  there  is  relatively  little  difference  in  trajectories 
as  (1  is  varied.  The  major  effect  of  increasing  3 is  the  faster  decay  of  shock 
strength  at  low  z and  slower  decay  at  large  z previously  shown  in  Figure  9. 

Tn  addition  there  appears  to  be  only  a very  weak  dependence  of  the  trajectory 
ijpuii  I us  Lh(!  curvis:  for  Y = l.l  fall  in  the  same  narrow  band  indicated  for 
i l.l*.  The  Y is  more  (‘vident  iti  the  z-y  plots  of  Figure  11  which  also 

iinlicato  the  e/Tect  ol'  a Cor  3 = 0.5  in  all  cases. 

TTie  apportionment  of  energy  between  kinetic  and  internal  modes  is  also  of 
interest.  These  two  contributions  can  easily  be  separated  in  the  conqjutation 
of  the  energy  integral,  J.  The  fraction  of  the  total  energy  associated  with 
kinetic  energy  is  plotted  in  Figure  12  versus  the  shock  strength  with  3 as  a 
parameter.  These  results  were  computed  with  a three  term  approximation  of  J. 
Note  that  the  kinetic  energy  fraction  increases  with  3 while  in  all  cases  most 
of  the  energy  soon  is  transformed  to  internal  energy  as  the  shock  decays  and  y 
Increases.  The  values  obtained  for  y = 0 are  in  agreement  with  computations 
made  by  Rogers^. 
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TABLE  1.  EXPANSION  COEFFICIENTS  FOR  CYLINDRICAL  WAVE 
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Flow  Field  Structure  Behind  Cylindrical  Blast  Wave  for  Constant  Energy  Input  Rate. 
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Decay  Index  and  Shock  Front  Position  from  Second  Order  Perturbation  Solution  (P) 
Compared  with  "Exact  Solution"  of  Goldstine  and  Von  Neumann  (G) 
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Sunmary  of  first  and  second  order  equations  required  for  blast  wave  computation 
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r<<««<<<«<<««<  blast  wave  *VAR  IABLE  ENERGY  »»>»»>»>>»»>»»>CBLST 


r.  PPOGRAH  COOED  BY  PROFESSOR  E.T,  PITKIN  CBLST 

MFCHANICAI  ENGINEERING  OEPARTNENT  CBLST 

IIMVEPSITY  CF  CONNECTICUT  CBLST 

STfPRS,  CONN.,  06268  CBLST 

«<«<<<<<««’<<<«_SECONJ)  ORDER  COMPUTATION  >»>»>»»»>»>»>»»>CPLST 
IMPL  IfIT  REAL*8  ( A-H,n-Z‘)  ' ' BLST 

. PEAL*9  MASS.lAMN.LAMl  ,lAM2tLAM3tlAM4,KN,KUt AHPI,LAMDA.LAMN2  BLST 

OIMENSICN  2A11 (200) «ZAP11( 200 ) «ZA12 1 200) t ZAP12 1200) * ZGl I ( 200 ) t BLST 
lZG12(200),ZA?I(200)f  ZA?2<200)«ZAP2H200),ZAP22I200)  »ZG2H200)«  BLST 

2ZG?2r?0C),C<  16),ZA(2bo)  tZAP(206)  .Z  Ml  200).  ZGV200)  * BZ  ( 100)  .0M|  loo)  BLST 
__  EQUIVALENCE  I C 1 1 ) . L AMN ) . IC ( 2 ) .KN ) , I Cl  3 ) .K1 ) . 1C  I ») . ALP) . I C ( 5) .GM)  , BLST 
1 (C  (6),GMQr.)  .ICI7)  «AK11  ) ,ICI6)  *AK21 ) • 1C  ) *GN)  . I Cl  10)*RLAM),  BLST 

_2  iri  ID.LAMl),  I C(  12)  .ALFA)  .1  Cl  13). GAMMA)  BLST 

CCNSTPUCT  INITIAL  VALUFS  FOR  INTEGRATION  BLST 

SAVE  * -.OOSDO  __  BLST 

A FnRMATI 3Fi6.«.I5)  BLST 

1Q_PEAD  (_5,<i,EN0*«)<i'))  ALFAjBETA^GAMMA.KFIELD  BLST 

L * ALFA  ♦ .00100  ' BLST 

LP  * 1*1  BLST 

ALP  = ALFA  ♦ l.DO  BLST 

GDI  < gamma  ♦ 1.00  _ BLST 

GM  « gamma  - 1.00  BLST 

GMCG  * GiV  C,  Ammj^  PL  ST 

AKli  » ?^00/GMnG  - GMaG/2.0b  BLST 

LAMN  * 2.U0*(ALP-fleT^A)/(2,O0*BETA|  . BLST 

LAmn?  * IAMN**2  BLST 

LABiPl  « LAMN  BLST 

KN  » ?.00*r.AMMA*(GM/GPL/ALP)**GAMMA/GPL  BLST 

Kl  s 0.00  BLST 

El  » 1.00  - IAMN/ALP/GAMMA  BLST 

F2  » l.nO  ♦ GMOGBL  AMNZALP  BLST 

F3  « I. 00  ♦ (7.00BGAMMA  - 1. 00) *LAMN/ GAMMA/ALP  BLST 

KPASf  » 1 BLST 

NV  » 12  BLST 

WPITEI6.A0)  BLST 

W0I^F(6.A1)  ALFA. BETA, gamma  BLST 

20  STEP  * SAVE  BLST 

KK  » OABSI .OlOO/STEP)  ♦ .100  BLST 

on  21  T«l.is  BLST 

21  21  I ) * O.UO  BLST 

MASS  * 1.00  BLST 

7 ( ) ) * 1 .00  BLST 

Zl?)  * GM/GPL/ALP  BLST 

Zl  S)  « 2.00/ALP/GPL  BLST 

J * IRl  BLST 
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7A12(J)  « O.DO  BLST  510 

2A?1(J)  * 0.00  BLST  520 

ZA?2( J)  * 0.00  BLST  530 

ZAPll (J)  » Z(5 ) BLST  540 

ZAPl2fJ)  * 0.00  BLST  550 

ZAP2l(J)  * 0.00  BLST  560 

ZAP22(J)  * 0.00  BLST  570 

ZG(J)  = KN/Z ( 2 )*«GAMMA  BLST  580 

ZGll(J)  * ZG(J)*(AKll  - GAMMA*Zf 5) /Z( 2) ) BLST  590 

ZG12(J)  = 0.00  BLST  600 

ZG2l(J)  = 0.00  BLST  610 

ZG22(J)  = 0.00  BLST  620 

25  W«)ITF(f,42}  ZA(  J),ZAP(  J)  tZC  31  .ZGC  J)  ,ZBCJ)  BLST  630 

IF(KPASS.E0.2)  GO  TQ  BLST  640 

WPTTE(6,43)  ZAIK  J)  .ZAPIK  J)  »ZA12(  J)  tZAPl2(J  )»  ZI8),  Z(9)  BLST  650 

IF(KPASS.EO.l)  GO  TO  50  BLST  660 

30  WRITE(6,44)  ZA? I ( J ) , Z AP21 ( J) ♦ZA22 < J ) t Z AP22CJ ) t Z C 14) , Z ( 15  ) BLST  670 

40  FORMAT!////)  BLST  680 

41  FOPWAT!*  ALFA  * • . F5  .3 « 3X  t *BETA  **tF  5.  3»  3Xt ’GAMMA  *»,F6.4,//)  BLST  690 

4?  F0PMAT(/,3X, *A0  =•  ,E10.4*3X,*A0PPIME  * • . E 10. 4, 3X t • JO  »».Ell.5,3X  BLST  700 

I, 'GO  *•  ,E1  1.5,  3X, ’MASS  «’,F10.3)  BLST  710 

43  F0PMAT(/,3X,’A11=’,E10.4,3X,’AIIPRIME«’,E10.4,3X,’A12«’,E10.4,3X  BLST  720 

l,’A12PRIME«’,E10.4,3X,’Jll  = ’,Ell.5,3X,’J12*’,E11.5,/)  BLST  730 

44  FnPMAT(/,3X,’A2l=’,E10.4,3X, ’A2lPRIMEx’,E10.4,3X, •A22«',E10.4,3X  BLST  740  ^ 

I, ’A22PRIME=’,E10.4,3X,’J2l*’,E11.5,3X,’J22*’,E11.5,/)  BLST  750 

45  FORMAT!/, » njO  = ’ , F 10 . 5, 2X , ’ EKO  * ’ , F 1 0. 5, 2X, ’EF RAC 0 x',F10.5)  BLST  760  j 

46  FORMAT!/,’  PLAM  =*  ,F I 0.5 ,2 X , • K2 I *’tFl0.5)  BLST  770 

47  FORMAT!/,’  OJll  = ’ , F 10 . 5, 2X, ’OJ 12  » • , F 10. 5,2 X, ’LAMl  x’,Fl0.5,2X,  BLST  780  ^ 

I’Jl  *’,F10.5,2X, ’EKl  *’,F10.5,2X,*KI  »’,F10.5)  BLST  790  ^ 

48  FORMAT!/,’  0J2I  *’ ,F 10.5,2X, ’Oj22  * ’ , FIO. 5,2X, ’ LAM2  x’,Fl0.5,2X,  BLST  800 

1’J2  =’,F10.5,2X, ’EK?  x’,Fl0.5)  BLST  810 

COMMENCE  INTEGRATICN  BLST  820 

50  00  100  1*1, KK  BLST  830 

100  fall  RUNKUT!Z,MASS,C,STEP,NV)  BLST  840 

J * j-i  BLST  850 

FMl  X MASS4*!LAMN/ALP)  BLST  860 

TF!KPASS.EO.?)  GO  TQ  98  BLST  870 

ZM!J|  X mass  BLST  880 

ZA! J ) * Z! 1 ) BLST  890 

ZAP! J)  * Z! 2)  BLST  900 

ZG!J)  * KN/! !Z!l)**L*Z!2) )**GAMMA*FM1 ) BLST  910 

ZAll!J)  * Z!4)  BLST  920 

ZA12!J)  * Z!6)  BLST  930 

ZAPIUJ)  * Z!5)  BLST  940 

ZAP12!J)  * Z!7)  BLST  950 

ZG11!J)  * ZG! J)*!AK11*FM1  - GAmma*! ALFA*Z !4) /Z ! 1 ) tZ  ! 5) /Z!  2) ) ) BLST  960 

ZG12!J)  * ZG! J )♦! ! 1,00-FM1|/LAMN-GAMMA*!ALFA*Z!6 )/Z ! I )tZ!7)/Z! 2) ) )BLST  970 

98  IF!XPASS.E0.1  ) GO  ’’O  99  BLST  980 

ZA21! J)  « Z!10)  BLST  990 

ZA??! J)  * Z! 1?  ) • BLSTIOOO 

ZAP?)  !J)  * Z!n  ) BLSTlOlO 

AAr??(j»  ■ Z<13)  PLST1020 

M.l  - ♦ t AM1*zr.)?!J)  BLST1030 

; A 1 - : « 4 I « ) A'tj  ♦ • (ft  ) «L ST  104  0 
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Zr,?l(J)  *(  ( (ZGl/ZGf  J)  )**2-RLAM**2»/2.D0  - GAMM  A*(  AL  FA*Z  ( 10) /Z  ( I ) 

1 4^  Z(lli/Z(2)  - (ALFA*(ZA1/Z(1))**2  ♦ ( ZAPl/ Z(  2 ) ) **2  ) /2 .00) 

2 ♦ AK21*FMl**2  - K l*RL AM*FM I )*ZG( J ) 

ZG22(J)  * ( (1  .D0-Ff^l**2)/2.00/LAMN  - GAMMAP<  AL  F A*Z<  1 2) /Z  ( I ) ♦ 

1 7(  13)/Z(2  ) ) )*ZG(  J) 

CUT  DOWN  STEP  SIZE  AS  MASS  APPROACHES  ZERO 
IFfJ.EO.lOl)  STEP  * STEO/IO.DO 
IF(J.EO.lOl)  GO  TO  101 
IF(J.EO.ll)  GP  TO  101 
1 F(J.GT  .?)  GO  TO  50 
101  Z3  « -Z(3) 

MRITF(6,42)  ZA(  J)*ZAP(  J)  tZ  3»ZG(  J)*ZM(  J) 

I FfKPASS.EQ.2  I GO  tq  104 
Z8  * -Z18) 

Z9  * -Z(P) 

WRITE  (6,43)  ZAIK  I),ZAP11  ( J),  ZA12(  J)  ,ZAP  12(  J)  ,Z8,  Z9 
IF  (KPASS.EO.l)  GO  TO  105 
10  4 714  * -Z(14) 

Z15  * -Z(  15) 

WRITF(6,44)  ZA21  (J) , Z AP2  1(  J ) ,ZA22(  J ),  ZAP22  ( J ) , Z 14,  Z 15 
105  IF(J,GT.2)  GO  TO  50 

COMPUTE  EXTRAPOLATED  VALUES  FOR  MASS  * 0 WITH  GO  * CONSTANT 
ENSTP  = -.0000100 
on  202  JK»  1,  2 

on  200  1*1, 90 

200  CALL  PUNKUTIZ, mass, C, ENSTP, NV) 

FMl  « MASS**(LAMN/ALP) 

GN  * KN/<(Z( 1)4*L*Z(2))*4GAMMA*FM1) 

73  * -Z(3) 

WRITE  (6,42)  Za),Z(2),Z3,GN,MASS 
IF(KPASS.E0.2)  GO  TO  199 
Z8  * -Z(8) 

79  * -7(9) 

WRITE (6, 43)  7(4),Z( 5) ,Z(6),Z(7),Z8,Z9 
IF(KPASS.tQ.l I GO  TO  201 
199  7 14  * -Z( 14) 

Z15  * -Z115)  . 

WRITE (6,44)  Z( 101 , Z(ll ), Z( 12),Z( 13) ,Z14,Z15 

..201  1F,UK.LQ^21.,GCLIQ.  ZQZ 

Z2S  « 7(2) 

7 4S  » 7(4J  _ . _ 

Z5S  * 7(5) 

Z65  * 2(6l) _ 

77S  » Z(7) 

ONS  * CN 

FNSTP  * ENSTP/10.00 
IF(XPA5S,E0.1)  GO  TO  20? 

7H1S  * 7(10) 

/ns  ’ /(ID 
/(.•••  - / ( I ; I 
/ I IS  - /(  D) 

’i>  • t itM  I NUt 

If  (KPAS5.E0.?)  GO  TO  206 
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Z(l) 

AN**L 

BLST1590 
RL ST  1600 

ANP  * 

Z(2) 

BLST1610 

FAJ  * 

AN  -ALP^ANPOMASS 

BLST1620 

GNO  * 

(GN/GNS)*GN 

BLST1630 

GN  « 

.7D0*GN  ♦ .3D0*GN0 

BLST1640 

GNKF 

« GNOlKN/GNlOPd  .OO/GAMMAl/GM 

BL ST  1650 

FJl  « 

GNKF*MASS**E1/El 

BLST1660 

FJ?  * 

GNKF*MASS*»E2/E2 

BLST1670 

FJ3  » 

GNKF*MASS**E 3/E3 

BLST1680 

ANO  « 

AN**LP  - ALP*(GM/GN)*FJ1 

BLST1690 

IF(AN0.LE. l.D-13)  AN0«1.0-13 

8L ST  1700 

ANO  « 

ANO**!!. DO/ALP) 

BLST1710 

‘ ANPO 

* (Z(2)/Z2S1*Z(2) 

BLST1720 

DJO  « 

GAMMA4'MASS4'FAJ4'*2/ALP/2.D0  ♦ FJl 

BLST1730 

A JO  * 

DJO  - Z(3) 

8LST1740 

EKO  * 

DJO  - FJl  - Z(10) 

BLST1750 

! EFPACO  « EKO/AJO 

BLST1760 

AMASS 

« O.DO 

BLST1770 

WRITE (6 *42)  ANO* ANPO* A JOtGNOt AMASS 

BL ST  1780 

WRITE(6«45)  OJOtEKOtEFPACO 

BLST1790 

All  « 

Z(4) 

BLST1800 

APll 

« Z(5) 

BLST1810 

A12  * 

Z(6) 

BLST1820 

AP12 

« Z(7I 

BLST1830 

Alio 

* (A11/Z4S)*A11 

BLST1840 

A 120 

» (A12/Z6S1*A12 

RLST1850 

APllO 

* (AP11/Z5S  )*AP11 

BLST1860 

AP120 

* (AP12/Z7S1*AP12 

BLST1870 

Gil  * 

GN*(AK11*FM1  - GAMMA4(ALFA«Z(4)/Z(l)«Z(5)/Z<2n) 

8LST1880 

G12  * 

GN*(  (1.00-FMl)/LAMN-GAMMA*<ALFA*Z<6l/Z<l)*Z(7)/Z<2n) 

BLST1890 

FKll 

» GAMMA*MASS*FAJ*(LAMPl*All/ALP  - MASS*AP11) 

BLST1900 

EK12 

» GAMMA4>MASS4FAJ*(LAMP1*A12/ALP  - MASS4AP12) 

BLST1910 

Dill 

» EKll  ♦ AK110FJ2  -GM*(ALFA*All/AN»APll/ANP)4FJl 

BLST1920 

0J12 

• EK12  ♦ (FJ1-FJ2)/LAMN  - GM*( ALFA*A 12/AN^AP12 /ANP) PFJl 

BLST193a 

EKll 

« EKll  - Zdl) 

8L ST 1940 

EK12 

* EK12  - Z<12) 

BLST1950 

! AZ8  « 

DJll  - Z(8I 

BLST1960 

A7M  • 

0J12  - Z(9) 

BLST 1970 

WRITE(6*43)  AllOf APllOf A120fAP120«AZ8fAZ9 

BLST1980 

ZM(1) 

* MASS 

BL ST 1990 

ZA(  1 1 

« AN 

8LST2000 

ZAP(l)  « ANP 

BLST2010 

1 ZG(l) 

* GN 

BLST2020 

1 ZAIKI)  « All 

BLST2030 

ZA12(  1)  * A12 

BLST 2040 

ZAPll(l)  « APll 

BLST2050 

ZAP12(1)  « AP12 

BLST2060 

ZGIKI)  * Gll 

BLST2070 

ZG12(l)  « G12 

BLST2080 

IFlKPASS.EO.n  GO  TO  210 

BLST2090 

206  WRITE(6f42)  ANOf ANPOt AJO f GNO t AMASS 

BLST2100 

A1  * 

All  * LAMIOAIZ 

BLST2110 

API  * 

APll  ♦ LAM1*AP12 

BLST2120 

T 


BBf  AVAILABIE  COPY 


A21  * ZilOl 
A22  * Z(12l 
AP21  « Z(ll) 

AP22  » Z(13l 
A210  « (A21/Z10S|4>A21 
A220  » <A22/Z12SI*A22 
AP210  » (AP21/ZllSI*AP2l 
AP220  * <AP22/Z13SI*AP22 
r,l  » Gil  ♦ LAM1*G12 

EK21  * GAMMA*MASS*(FAJ*(  A21*aAHPl*LAMNI  - ALP4'HASSOAP21  ♦ LAMl 
1*A1I  ♦ <LAMP1*AI  - ALP*m'ASS*AP1  l♦♦2>2.^00^/ALP 
0J21  « EK21  *■  AK21*FJ3  - K1*RLAM*FJ2  ♦ F Jl4>(  (G1 /GN)  «>(  ALFAOAl  / AN  ♦ 
1 APl/ANPI  ((G1/GNI««2  - RLAM**2  ♦ G AMHA4>  ( ALF  A^f  Al/AN  1 4'4'2  ♦ <AP1 
2/ANP - GM*<ALFA*A21/AN*AP21/ANP)  ♦ 

3 -l.D0l*Al'/AN>2.00'*APl/ANP)) 


EK22 

r)J22 

EK21 

EK22 

AZ14 

AZ15 


EK21 
EK22 
0J21 
OJ22 
WRITE(6,44I 
ZA21 (II 
ZA22( 1) 
ZAP21  (1  I 


GAMMA*HASS*FAJ*(A22*  (LAMPl+LAMNI  - ALP*M ASS*AP 221 /ALP 
EK22  *■  (FJ1-FJ3I/LAMN/2.00  - GH*  ( ALFA*A22/AN*AP22/ANP  I ♦F  J 1 

- Z(8I 

- Z(9I 

- Z(14) 

- Z(15I 

A210.AP210f  A220tAP220fAZ14tAZ15 
A2l 
A22 
AP21 


/AP22<1I  » AP22 

/G21(ll  « (((G1/GN|**2  - RLAM**2l/2,()0  - GAMMA*  ( ALF  A*  Z ( 10) /Z  ( 1 1 

1 * Z(11)/Z(2I  “ (ALFA*(  Al/Z(l)l**2  ♦ ( AP 1/Z ( 2)1 **2 I /2.00) 

2 AK21*FM1**2  - K 1*RL AM*FM1 ) ♦GN 

ZG22(1)  * ( ( l,00-FMl**2)/2.00/LAMN  - GAMMA*( ALFA*Z( 12) /Zl 1)  ♦ 

1 Z< 13)/Z(2) ) )*GN 
IF(KPASS.E0.2)  GO  TO  220 

COMPUTE  LAMBDAS  AND  ENERGY  INTEGRAL  CONTI  BUT  IONS 
210  ONEGO  * 2«00/(LAMN4-2.00) 

OMEGl  ■ (2*00*LAMN«-2.DO)/(3*00*LAMN+2.00) 

LAMl  * (1,00/GM/8LP-AZ8)/(AZ9+AJO*(BETA*OMEG1/OMEGO-ALP)/LAMN2) 
RLAM  * LAMl/LAMN 
K1  a AKll  - RLAM 

A-  21  * AK11*RLAM  - 2.00/GM0G/GM  - GM0G**2/8.00  - RLAM**2/2.00 

Aj.  * AZ8  ♦ LAM1*AZ9 

EKl  EKll  *■  LAM1*EK12 

WRITE(6t47)  OJllrOJUrLAMltAJlfEKltKl 

WRITE(6f46)  RLAM«AK21 

NV  « 15 

KPASS  * 2 

GO  TO  20 

220  BZl  . -LAM1/LAMN2 

0MEG2  « (4, D0*LAMN>2. 001/(5. D0*LAMN*2.00) 

SA  « (8ETA*0MEG2/0MEG0  - ALP)/2.nO 
SB  » -AaO«A/L\MN2  _ 

SA  «AJ0*8Z1**2*( SA*LAMP1  ♦ BETA* ( BETA>1. 00 )*( OMEGl/ OMEGO) **2 
_ 1/2^00  ♦ ALP*(  (ALFA»2.D0)/2.00  - BFTA*0MEG_1/0MEG0| ) 

LAM2  « (SA  - AZ14)/(AZ15  - SB) 

LAM4  «-(  l.OO-'BETA/ALP)  ♦((  ALFA+2.00) /4.00  - LAM2  - 2.00*LAM1  - 


BLST2130 
BLST2140 
BLST2150 
BLST2160 
BLST2170 
BLST2180 
BL ST2190 
BLST2200 
BLST2210 
BLST2220 
BLST2230 
BL ST  2240 
BLST2250 
BLST2260 
BLST2270 
BLST2280 
BLST2290 
RLST2300 
BLST2310 
BLST2320 
BLST2330 
BLST2340 
BLST2350 
BLST2360 
BLST2370 
RLST2380 
BLST2390 
BLST2400 
BLST2410 
BLST2420 
BLST2430 
BLST2440 
BLST2450 
BL  ST  2460 
BLST2470 
BL ST 2480 
BLST2490 
BLST2500 
BLST2510 
BLST2520 
BLST2530 
BLST2540 
8LST2550 
BLST2560 
RLST2570 
RLST2580 
BLST2590 
BLST2600 
BLST2610 
BLST2620 
BL  ST  26  30 
8LST2640 
8LST2650 
BL ST 2660 
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Bffl'AVAIUBlE  COPY 


I 3.00>*LAHN4 
LAM  3 » -LAMN  - 
AJ2 
-LK?. 


iAJlJL^JLAMl, -.i.AM4_ 


AZ14 

^.K2I 


LAM2*AZ15 

LAM2*EK22 


BLST2670 
RLST2680 
BLST2690 
BLST2700 
BLST2710 
BLST2720 
BLST2730 
BLST2740 
BLSf2750 

BLST2760 

-(  LAM3*  ILAMPl*LAMNl*LAMl*BZI2»+LAMPm.  AM2*B2  II/3.00/LAMN2  BLST2770 
~l LAM4»( l.QQ*3.DO*LAMN)*LAMl*BZI3)»ILAMPl»LAMNI*LAM2*BZt 2)BLST  2780 

8LST2796 
BLST2800 
BLSt28l'6 
BLST2820 
BLST2830 


WRITE  16 *48)  OJ21*OJ22*LAH2*AJ2«EK2 

COMPUTE  SHCCK  TRAJECTORY ____  _ 

WRITEI6,40) 

BZQ  » U(2.DQ»BETA)/(ALE;^>3.DQ))**^TA/AJQ)»*J_UD0/  < AlP-BfTA  ) ) 

RZIl)  * BZl 

B7C2)  X -<LAM2  » LAMP  1 AM1*BZ  1 1 /Z>Q0/LAMflt2. 

BZI3) 

- BZ(41 


1+L AMP1*LAM3*BZI )/4.D0/LAMN2 

■ QQ  221  K«g.*lQO 

AK  s K 

BZ(K)  = -I  ( l.DQ»<AK-4.(X))*LAMN)*8ZIK-4)»LAM»  ♦ ( l.OOMAK-3.00) 


l*LAMN)*BZ(K-3)*LAM3«-(  1.00«-l  AK-2.  00) ♦LAMN)  *B2(K-2 ) AM2 


2-  * (1.00»(AK-I.00 J ♦LAMN )^BZ(K-1) ♦LAM1)/LAMN2/ AK 

221  6m(K)  = 1. 00/(1. DO*. 500/IAK4-1.D0/LAMN)) 

_ QM(1)  « OMEGl 

0M(2)  * 0MEG2 

0M(3)  « 1.0Q/(  1.00». 500/(3. 00*1. OO/LAMN)) _ 

0M(4)  * 1.00/( 1 .00*. 500/(4. 00+1. 00/LAMN) ) 

a.DjQ...  . - - 

00  225  I«l,20 

SDMZ  « Q.QQ 

SUMO  * 0.00 

00  222  K»l,100 

TEPMZ  * BZ(K)^Y^^K 
TEPMC  * 0M(K)4TERMZ 

IE(TERMO.LT.l. 0-5. ANO.TERMZ.LT. 1.0-5)  GO  TO  223 
SUMZ  * SUMZ  ♦ TERMZ 

222  SUMO  * SUMO  ♦ TERMO 

223  EIRST  * BZO*Y^^( 1.00/LAMN) 

BZS  « (1.00  ♦ SUMZ)^FIRST 

OMEGA  * (OMEGO  ♦ SUMO  ) ♦F IR  ST^OSORT  ( Y ) 

LAMOA  » LAMN  ♦ LAMl^Y  ♦ LAM2^Y^^2  *■  LAM3^Y^^3  ♦ LAM4^Y^^4 
AJ  * AJO  ♦ AJl^Y  ♦ AJ2^Y^^2 
EK  ■ EKO  ♦ EKl^Y  ♦ EK2^Y^^2 
EFPAC  » EK/AJ 

WRITE (6 *440)  Y* OMEGA  * BZS *AJ * EFR AC *LAMOA 
225  Y » Y ♦ .0500 

WRITE (6 *442)  ALFA*8ETA*AJ0* AJl *AJ2*LAMN*LAMl *LAM2*LAM3*L AM4* 
1 AZ8*AZ9*AZ14*AZ15*GAMMA*SAVE 
:OMPUTE  FLOW  FIELD 

IF (KFIELO.EO.O)  GC  TO  10 
00  350  II»1*9 
Y » o.ino^( II-l) 

WRITE(6*40) 

LAMOA  * LAMN  ♦ LAMI^Y  ♦ LAM2^Y^^2  ♦ LAM3^Y^^3  *■  LAM4^Y^^4 
WRITE(6*445)  Y*LAMDA 
WRITE (6*446) 

DO  350  KK*1*110 
IF(KK.EO.l)  GO  TO  302 


BLST2840 
BLST2850 
BLST2860 
BLST287b 
BLST2B80 
BL  ST  2890^ 
BLST2900 
BLST2910 
BLST2920 
BLST2930 
BLST2940 
BLST2950 
BLST2960 
BLST2970 
6LST2980 
BLST2990 
BLST3000 
BLST3010 
BLST3020 
BLST3030 
BLST3040 
BLST3050 
BLST3060 
BL ST3070 
BLST30P0 
BLST3090 
BLST3100 
BLST3110 
BLST3120 
BLST3130 
BLST3140 
BLST3150 
BLST3160 
BLST3170 
BLST3180 
BLST3190 
BLST3200 


BEST  AVAIIME  COPY 


IF(HCO( (KK-l)»5).NE.O. ANO.KK.lt. 1001  GO  TO  350 
IF  (KK.LT.92I  GO  TO  302 
IFIKK.GT.IOOI  go  to  301 
JJ  « lOl  - (KK-9ll*lO 
GO  TO  303 
301  JJ  » III  - KK 
GO  TO  303 
30?  JJ  » 192  - KK 

303  ZAl  * ZAIKJJ)  LAM1*ZA12(  JJI 
ZA2  « ZA21(JJ)  * LAM2*ZA22(JJI 
ZAPl  * ZAPIKJJ)  4^  LAM1*ZAP12(JJI 
ZAP2  » ZAP21(JJ)  ♦ LAM2*ZAP22( JJ) 

ZAT  * ZACJJI  ♦ ZAl*Y  ♦ ZA2*V**2 
ZAPT  » ZAP(JJ)  ♦ ZAPl*Y  ♦ 7AP2*V**2 
V » ALP*ZAT**L*ZAPT 

VEL  * ZAT  - ALP*ZM(JJ)*ZAPT  ♦ LAMOA*Y*(  ZAU2.D0*ZA2*Y) 

ZGT  < ZG(JJ)  * (ZGll(JJ)  4^  LAH1*ZG12(  JJn*Y  4^(  ZG2  1 ( J J 1 4^LAH2« 
I ZG22(JJ)I*Y**2 


BLST3210 

BLST3220 

BLST3230 

BLST3290 

BLST3250 

BLST3260 

8LST3270 

BLST3280 

BLST3290 

BLST3300 

BLST3310 

BLST3320 

BLST3330 

BLST3390 

8LST3350 

BLST3360 

BLST3370 

8LST3380 


IFCKK.GT.ll  GO  TO  305 
ZGTF  « ZGT 
VF  * V 
VEIF  » VEL 
30  5 RHO  » VF/V 

ZGT  » ZGT/ ZGTF 
VELR  » VEL/VELF 

WRI TE(6,44B)  ZM ( J J | , Z A 7,ZGT , RMOt VEL R 
350  CONTINUE 
GO  TO  10 
440  FORMAT!  • 
l.‘  KE/J 


Y **,F6.3,*  OMEGA  »*,Ell.5,* 

'*,F10.5,*  LAMOA  »»,F10.5) 

442  FORMAT!/ 

lFl0.5,2Xi • JI**,F10.5, 2X, *J2»* ,F10.5,/,'  LA MO*' tFlO. 5,2X , »LAMl*' , 


SUMMARY'  ,/,'  ALFA*'  tFl0.5»2X,  'BETA»'tFl0.5*  2X»  'JO*'  t 


2F10.5,2X,' LAH2*'  ,F  10. 5 , 2X,  ' L AM3* ' t F 10  .5,  2X , ' LAM4*  ' , F 10.  5,  /,  ' Jll* 


BLST3390 
BLST3400 
BL ST 3410 
BLST3420 
BLST3430 
BLST3440 
BLST3450 
BLST3460 
BLST3470 
BLST3480 

Z •'♦Ell. 5,'  J ••♦FIO.5  BLST3490 

8LST3500 
BLST3510 
BLST3520 
'BLST3530 


3,FI0.5,2X, 'J12*' fFl0.5t2X,' J?l*' .FIO. 5 ♦ 2X ♦ ' J22*' tFlO.St/f'  GAMMA*' RLST3540 


4,Fl0.5f  2X,  'STEP*  '♦  FIO. 51  BLST3550 

445  FORMAT!/,'  Y *' ,F5.2 tZX, ' LAMOA  »',F10.5)  BLST3560 

446  FORMAT! 4X, 1 5X, 'MASS ', 13X , 'RADIUS '♦ 1 1 X, ' PRE SSURE *♦ 12X  ♦•DENSITY'  BLST3570 

l,llXf'VELOCITY't//l\. _ BLST3580 

448  F0RMATI4X, 5F19.5I  BLSt3590 

999  STOP  BLST3600 

END  BLST3610 


BESI'AVAIUBLE  COPY 


SUBROUTINE  FCT( Z ,ZP,C , N| 

IMPLICIT  REAL*8  (A-H,0-Z) 

PEAL*8  LAMN,KN«M,LAM1,K1,LAMP1 
DIMENSION  Z( 15) «ZP( 15) •C(16) 

C AND  Z TRANSFERS 

IF(M.NE.l.DO)  GO  TO  10 
LAMN  « C(l) 

KN  » C(2) 

K 1 « C<  3) 

ALP  « C(4) 

GAMM  » C<  5) 

GMCG  * C<6) 

AKll  « C(7I 
AK21  « C<8) 

RLAM  » C(IO) 

LAMl  * C(ll) 

ALF  « C(12l 
GAM  » C(13) 

LAMPl  > LAMN  * l.DO 
L * ALF  ♦ .0100 
10  AN  > Z( 1) 

ANP  « Z(2) 

COMPUTE  COMMON  TERMS 
ANL  > AN**L 
FMl  « M*«(LAMN/ALP) 

GN  « KN/(ANL*ANP)**GAM/FM1 

FAJ  « AN  - ALP*M*ANP 

AA  »_  ALP*GAM*(ANL*GN/ANP_-_ALP*m**2) 

COMPUTE  ZERO  ORDER  DERIVATIVES 

B « ALP*GAM*(2.00*ALF  ♦ LAMN)*MZ2.00 

CO  » GAM*(  LAMN/2.00  ♦ GN *AL F*ALP*ANL*ANP/AN**2 ) 

ZPai  « ANP 

ZP(2)  « (8*ANP  - Ca*AN  - GNRANL*LAMN/M ) /AA 
EKO  ■ GAM*FAJ**2/ALP/2.00 
ZP(3)  « EKO  ♦ GN*ANL*ANP/GAMM 
COMPUTE  FIRST  ORDER  DERIVATIVES 
All  - Z(A) 

APll  ■ 7(5) 

A12  « Z(6) 

AP12  - Z(7) 

ANPP  « ZP(2) 

GNP  » -GN*(LAMN/ALP/M  ♦ GAM*{ALF*ANP/AN+ANPP/ANP) ) 

BB  ■ ALP*GAM*(  ANL*(GN*ANPP/ANP**2  - GNP/ANP  - ALF*GN/AN) 
1 -2.00*LAMN*M)  ♦ B 

CC  ■ CC  ♦ GAM*LAMN**2/2.D0  - GAMM*ALF ♦ALPPANLPGNP/AN 
DO  « AK11*FM1*ANL*(GNP*ALP  ♦ GN*LAMN/M) 

EE  * (ANL*ALF*GNP  - DO/ AK 1 1 ) /L AMN  - FAJ*GAM/2.D0 
ZP(4)  « APll 

7P(5)  « (8B*APll+CC*All+00)/AA 
ZP(6)  > AP12 

7P(7)«  (BB*AP12*CC*A12*EE)/AA 
IF( K1 .NE.O.DO)  GO  TO  100 

EKll  « GAM4FAJ«( A11«LAMP1-ALP*M*AP11 )/ALP 
FK12  « GAM*FAJ*(A12*LAMPI-ALP*M*AP12)/ALP 


FCT  10 
FCT  20 
FCT  30 
FCT  40 
FCT  50 
FCT  60 
FCT  70 
FCT  80 
FCT  90 
FCT  100 
FCT  110 
FCT  120 
FCT  130 
FCT  140 
FCT  150 
FCT  160 
FCT  170 
FCT  180 
FCT  190 
FCT  200 
FCT  210 
FCT  220 
FCT  230 
FCT  240 
FCT  250 
FCT  260 
FCT  270 
FCT  2M 
FCT  290 
FCT  300 
FCT  310 
FCT  320 
FCT  330 
FCT  340 
FCT  350 
FCT  360 
FCT  370 
FCT  380 
FCT  390 
FCT  400 
FCT  410 
FCT  420 
FCT  430 
FCT  440 
FCT  450 
FCT  460 
FCT  470 
FCT  480 
FCT  490 
FCT  500 
FCT  510 
FCT  520 
FCT  530 
FCT  540 


i 

i 

j 
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BEST'AVAIWBIE  COPY 


ZP  (81  » EKll  ♦ GN*4NL*(AKll*ANP*FMl/t;‘AMM-APll-ALF*ANP*All/ANI 
ZP(9)  » EK12  ♦ GN*ANL  ♦(ANP*(  l.OO-FHD/GAHMZLAMN  - AP12 
1 - ALFA'ANP*A12/ANI 

2P11QJ  » EKO 

ZPdl  I » EKll 

?Pi-I2l  » EK12  ..  ...  

RETURN 

COMPUTE  SECOND  ORDER  DERIVATIVES 

100  A1  » All  ♦ LAM1*A12 

A.Pl  * APll  ♦ I.AM.I*AP12  _ _ 

APPl  » ZP(5)  ♦ LAM1*ZP(7I 

Gl  - GN*(RLAH  - GAM* ( ALF*A  1/ AN  _♦ . APl/ANPl  ♦ K I *FM I I 

GPl  » Gl*GNP/GN  ♦ GN*(KI*LAMN*FM1/H/ALP  - GAM*( ALF* ( AP l/AN-A  I 
l*Af!lP7AN**2|  ♦ ( APPJ/ANP-APl*ANPP/ANP**2)  n 
A 21  * Z( 101 


A22  » 2(121 

_A..p.I^2-JL  in3J . . 

BRB  ■ BB  - 2.00*ALP*GAM*M*LAMN 

CCC_  i CjC  ♦ 5A^EAMNJLU«D0*2.500*LAMISLI 

000  « AK21*LAHN*FMI**2*(GNP/IaMN+2’.D0*GN/ALp7m)  - Kl*LAMl*FMl* 
1 i GNP / LAMN»GN/ALP/M1  ♦ GNP *(GAM*(ALF ♦(Al/ANl** 2tJAPl/ANP)**2)  - 
2(  Gl/GNI**2-RLAM**2J/2.00+GPl*(Gl/GN+ALF*Al/AN»+GAM*GN*( AL F*Al* 

PPI-AP1**2*ANPP)/ANP**3) 


□ 00  * 000*ALP*ANL  ♦ GAM*L AM  I ♦( Al*( 1 , DO*4.00*LAMN » - 3.00*ALP*M* 

1 AP11/2.DQ  » ALP*ALF*(ALF-l«00l»Al**2»GNP/2.00  _ _ „ 

EFE  * .500*ALP*(  ANL*GNP«(1.00-FM1**2>/LAMN  ♦ 'GAM*M*ANP  ) - 

1 ANL*GN*FM1**2/M  - GAM*AN/2.00  

ZPdOJ  « AP21 

ZPdll  « (ReB»AP2 1»CCC*A2U00D> / AA 

ZP( 12  I * AP22 

ZP(13I  « ( BBB*AP22»CCC*A22*EEEI/AA  

G21  *(( (G1/GNI**2-RLAM**2|/2.D0  ♦ AK21*FM1**2  - RLAM*K1*FMI 
1 - GAP*(ALF*(A21/AN-( Al/AN ^♦*2/2.00^  ♦ AP21/ANP  - 1 AP I / ANP) ♦ ♦ 2 
2/2.00))*GN 

G22  ^ GN*( (1.DQ-FM1**2I/LAMN/2,00  - GAM*( ALE*A22/AN» AP 2 2/ ANP 1 1 
EK21  * GAMOlFAJ^lLAMl^Al  ♦ (L  AMPl+LAMN  )*A?l  - ALP*MA>AP2.l| 

1 ♦ (LAMP1*A1  - ALP*M»APll**2/2.D0>/ALP 

ZPdA)  * EK21  ♦'aNL*ANP*(  (AL  F*A2l/AN  * AP’2l/ANP  ♦ ( API/ANP+ ( ALF 
1-1.QQ1*A1/AN/2.001*ALF*A1/AN1*GN»G1*( AP 1/ ANP»ALF*Al / AN ) »G21 ) /GAMM 
EK22  » GAM*(F AJ*( (LAMP l+LAMN ^♦A22  - AL P*M*AP22 1) /AL P 
ZP(15I  * EK22  ♦ ANL*ANP*(GNO(ALF*A22/AN«-AP22/ANP)  ♦ G22)/GAMM 
ZP(8|  * EK21 
ZP(9I  * EK22 
RETURN 
END 


FCT 

550 

FCT 

560 

FCT 

570 

FCT 

560 

FCT 

590 

FCT 

600 

FCT 

610 

FCT 

620 

FCT 

630 

FCT 

64  0 

FCT 

650 

FCT 

660 

FCT 

670 

FCT 

680 

FCT 

690 

FCT 

700 

FCT 

710 

FCT 

72  0 

FCT 

730 

FCT 

740 

FCT 

750 

FCT 

760 

FCT 

770 

FCT 

780 

FCT 

790' 

FCT 

800 

FCT 

810 

FCT 

820 

FCT 

B30 

FCT 

840 

FCT 

850 

FCT 

B60 

FCT 

870 

FCT 

880 

FCT 

890 

FCT 

900 

FCT 

910 

FCT 

920 

FCT 

930 

FCT 

940 

FCT 

950 

FCT 

960 

FCT 

970 

FCT 

980 

FCT 

990 
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I 


BESF'AVAIIABIE  COPY 


SUBRCLTINE  RUNKUTf Z,X«C«H,N) 

RNKT 

10 

CONSTANT  STEP*  FOURTH  ORDER  RUNGE-KUTTA  INTEGRATION  ALGORITHM 

RNKT 

20 

implicit  REAL*8  CA-H,0-Z) 

RNKT 

30 

DIMENSION  Zf ISItCf 16) ,YI 15) ,RK1 ( 15 ) , RK2 f 15 ) t RK3 f 15) ,RK4fl5) 

RNKT 

40 

H2  « H/2.00 

RNKT 

50 

call  FCTfZ«RKl,C«X) 

RNKT 

60 

DO  10  J«1«N 

RNKT 

70 

10 

Yl J)  « ZIJ ) ♦ H2*RKlf J) 

RNKT 

80 

X * X*  H2 

RNKT 

90 

call  FCTfY«RK2,C« X) 

RNKT 

100 

DO  11  J«1«N 

RNKT 

110 

ll 

Y(J)  « ZfJ)  *■  H2*RK2fJ) 

RNKT 

120 

CALL  FCT(Y«RK3,C«X) 

RNKT 

130 

DO  12  J«ltN 

RNKT 

140 

12 

YU)  » ZfJ)  ♦ H *RK3fJ) 

RNKT 

150 

X « X^  H2 

RNKT 

160 

CALL  FCTlYtRKAtCtX) 

RNKT 

170 

DO  13  J«1«N 

RNKT 

ISO 

13 

ZU)  « ZfJ)  * H*fRKlf  J)'»2.00*fRK2f  J)4RK3f  J))»RK4t  J))/6.00 

RNKT 

190 

RETURN 

RNKT 

200 

END 

RNKT 

210 
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